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ABSTRACT 

Within the NUclei of GAIaxies project we have obtained IRAM PdBI and 30m '^CO(l-O) and '-C0(2-l) observations of the spiral 
galaxy NGC6951. Previous work shows that there is indirect evidence of gas inflow from 3kpc down to small radii: a large-scale 
stellar bar, a prominent starburst ring (rR;580pc) and a LINER/Seyfert 2 nucleus. In this paper we study the gas kinematics as traced 
by the CO line emission in detail. We quantify the influence of the large-scale stellar bar by constructing an analytical model of the 
evolution of gas particles in a barred potential. From this model gravitational torques and mass accumulation rates are computed. We 
compare our model-based gravitational torque results with previous observationally-based ones. The model also shows that the large- 
scale stellar bar is indeed the dominant force for driving the gas inward, to the starburst ring. Inside the ring itself a nuclear stellar oval 
might play an important role. Detailed analysis of the CO gas kinematics there shows that emission arises from two co-spatial, but 
kinematically distinct components at several locations. The main emission component can always be related to the overall bar-driven 
gas kinematics. The second component exhibits velocities that are larger than expected for gas on stable orbits, has a molecular gas 
mass of 1.8x10*Mq, is very likely connected to the nuclear stellar oval, and is consistent with inflowing motion towards the very 
center. This may form the last link in the chain of gas inflow towards the active galactic nucleus in NGC 695 1 . 

Key words. Galaxies: individual: NGC 6951 - Galaxies: nuclei - Galaxies: ISM - Galaxies: Seyfert 



^ 1. Introduction (ISM) are a hydrodynamical mechanism for losing angular mo- 

mentum. When present, gravitational torques are more efficient 
Gas accretion onto supermassive black holes (SMBHs) is be- g Combes . 2004 ) 
^-H lieved to be the cause of nuclear activity in galaxies. SMBHs 
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are accepted as a common co mponent in most galaxie s with a Observations of the inner regions of galaxies are neces- 
significant massive bulge (e.g. [Ferrarese & Ford||2005| and ref- sary to understand how the interplay of mechanisms resuks in 
erences therein). However, only about half (43%) of the local gas transportation all the way down to the nucleus. The dy- 



. ^ galaxies host active galactic nu clei (AGN) in t he Seyfert, LINER namic timescales in the central region of galaxies are short. 
5^ or transition-object categories (lHoetalJll997|). This discrepancy Therefore, the gas distribution in these regions needs to be 
9 ; between the presence of SMBHs and nuclear activity must there- mapped with high angular resolution. The NUclei of GAIaxies 
fore be sought in the possibility and efficiency of gas transporta- (NUGA) project (Garcia-Burillo et al. 2003a) has been obtain- 
tion to the central regions. Gas transport toward the centers of ing high-resolution (0.5" - 1") detailed mapping of the molec- 
galaxies can only happen when the gas is able to lose its angular ular gas kinematics in 12 nearby (D = 4 - 40 Mpc) low lu- 
momentum. Two categories of dynamical mechanisms can cause minosity active galactic nuclei (LLAGN) with the IRAM PdBI 
inflow. Gravitational mechanisms such as galaxy-galaxy inter- (piateau de Bure Interferometer) and 30m telescope. This sam- 
actions on lai-ge scales or non-axisymmetries (i.e. spiral den- pig spans the whole sequence of nuclear activity types. In the 
sity waves, large-scale stellar bars or nuclear ovals) within the central kiloparsec most of the gas is in the molecular phase, 
galaxy potential exert gravitational torques. Viscosity torques making CO lines optimal tracers of the gas dynamics. The spa- 
and shocks caused by turbulence in the interstellar medium tial resolution (< 100 pc) of fliis survey allows one to ob- 
serve the gaseous distribution over an impressive spatial range. 



* Based on observations carried out with the IRAM Plateau de Bure This has already led to the identification of a wide range of 
Interferometer and 30m telescope. IRAM is supported by INSU/CNRS morphologies in the nuclear regi ons of these galaxie s, includ - 
(France), MPG (Germany) and IGN (Spain). ing lopsided disks (NGC 4826: iGarcfa-BuriUo et all (l2003bh . 
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NGC 3 718: iKrit^s etalJ (l2005h . NGC 5953: [ Casasola et alJ 
(l2010h '). bars and spirals (NGC 4569: 'Boone et al." (2001), 
NGC 2 782: fHuntetal. (2008), NGC 6574: Lindt-Krieg et af, 
(l2008h . NGC 4579 Garcia-Burillo et all (120091)) and "rings 
(NGC 7217: Combes et al. 
l2008h . NGC 1961: Combes 



( 2004 ), NG C 3 147: ICasasolaltaLl 
et alMQOi)- 

Large-scale stellar bars are believed to be efficient in 
driving gas towar ds the i nner Lindblad r e sonanc e (ILR) (i.e. 
Regan & Teube3 12004 ISakamoto et aP Il999t IShethetalJ 
2005i) . There they induce the formation of spiral or 



ring structures (fPrendergast 1983; Athanassoula 1992a.b; 
Englmaier & Shlosman 2000; Macieiewski 2004; Combes 
2004). M artini & Pogge (1999 ) and RegaiL& Mulchaey ( 1999) 
have proposed that spiral shocks induced by the large-scale 
stellar bar can generate further gas inflow across the ILR. 
Numerical simulations have shown that spiral structure can 
extend across th e ILR into the nuclear region, if the sound speed 
is high enough dEngknaier & Shlosman 20 00) or if the ve locity 
dispersion in the ISM is large enough ?Macieiewskill2004l) . 

Alternatively, other mechanisms on smaller scales could take 
over gas transpor tation towards the nucleus, i.e. the bars-within- 
bars scenario (Shl osman et al.llT989h . Viscosity torques can only 
become the do minant mechanis m for inflow in the innermost re- 
gions (<200pc: ICombesI (l2004 ). 

In this paper we study the nearby spiral galaxy NGC 695 1 . 
Its distance of 24.1 Mpc allows for high spatial resolution (1" = 
1 17pc). There are several arguments that gas inflow is currently 
occurring in this galaxy. A large-scale stellar bar has been de- 
tected in the near-infrared (Buta et al. 2003; B lock et al.1l2004 
with a bar radius of 26" (a;3.0kpc) and a position angle of 84° 
(iMulchaev et al.l [1997) . This galaxy also shows a pronounced 
starburst (SB) ring at 5" (=580 pc) radius in Ho - (iPerez et al.l 



120001: iRozas et al. 



sion. iHaan et al 



l2002h and radio (ISaikia et all l2002h emi~ 
2009h have shown that the inner region of 



NGC 695 1 is HI depleted, implying the ISM must be dominated 
by its molecular phase. But '^CO(l-O) and '^C O(2-l) emis- 
sion associated with the SB rin g has been found (iKohno et alJ 
ll 999"; 'Garci'a-Burilloetal]|2005t), as well as HCN(l-O) emission 
(Krips et al. 2007). NGC 6951 has also been classified as 'a high 
excitation LI NER and a possibly high nitrogen abundant Seyfert 
T galax y (Perez et al l |200Q). In the central I" of NGC 6951 
iKrips et al. (2007) and| ' 



Garcfa-BuriUo et al.l (l2005h have detected 



HCN(l-O) and '^CO(2-l) emission, respectively, indicating that 
further inflow of gas beyond the SB ring must be or has been 
occuring recently. 

We present '^CO(l-O) and '^CO(2-l) observations made 
with the IRAM PdBI and 30 m telescope. Previous papers by 
iGarcfa-Burillo et al.l(l2005l) and lHaan et all (1200 9) presented the 
PdBI-only data. Here the '^CO(l-O) and ^^00(2-1) data cubes 
have been combined with the 30m observations. The addition 
of the 30m observations provides the full CO emission present, 
sampled on all scales in the inner 3 kpc of NGC 695 1 . 

This paper has three aims. The first is to investigate how 
the addition of the 30m data ch anges the maps and t he gra vita- 
tional torque resu lts derived by iGarcfa-Burillo et all (l2005h and 
iHaan et al] (|2009|) and to perform a more detailed investigation 
of the kinematics of the molecular CO gas than has been done 
before. The second is to quantify the influence of the large-scale 
stellar bar within the inner 3 kpc using a parametric kinematic 
model. The third is to study the impact of the nuclear stellar oval 
on the gas flow inside the circumnuclear gas ring. 

In §2 we present the observations from the IRAM PdBI and 
30m telescope and their reduction. §3 contains a discussion of 
the changes due to the addition of the 30m data and a presenta- 



tion of the spatial and kinematic properties of the CO emitting 
gas as seen in the PdBI+30m data cubes. In §4 we detail the 
large-scale stellar bar model, with the results that follow from 
the model being presented in §5. Finally, in §6 we compare our 
result with previous gravitational torque studies of NGC 6951 
and discuss observational evidence for inflow to the nucleus. We 
conclude in §7. 



2. CO observations 

2.1. IRAM PdBI observations 

The IRAM PdBI observations in ABCD configuration were 
carried out between June 2001 and March 2003 using the 6- 
antenna array in dual-frequency mode. Only the D-configuration 
observations were executed with 5 antennas. The correlator 
was centered at 114.726 GHz and 229.448 GHz at 3mm and 
1mm, respectively, corresponding to a heliocentric velocity of 
1425 km s '. The bandwidth covered changed between the CD 
and AB configurations, although all four configurations cov- 
ered at least the central +200 km s"^ The flux calibration used 
CRL618, 0932H-392, 3C273 and/or 3C345. Bandpass correction 
was derived from observations of a strong quasar at the begin- 
ning of the track. 1928+738 and 2010+723 served as the phase 
calibrators, allowing for correction of atmospheric effects. Phase 
corrections derived for the 3mm receiver were applied to the 
1mm band resulting in a better phase correction at 1mm. The 
data were reduced using standard routines in the GILDAS soft- 
ware packag43. 

For both data sets CLEANed data cubes were produced with 
natural and robusj^ weighting. Cleaning was done down to the 
2cr noise level within a fixed polygonal area that was defined 
based on the zeroth moment map for all channels with line emis- 
sion. The r.m.s. noise for these data cubes is listed in Table [1] 

2.2. IRAM 30m observation 

30 m observations of the central 132" by 66" were obtained 
on December 24 and 25, 1997. The 3mm and 1mm receivers 
were tuned to 114.730GHz and 229.460GHz. A bandwidth of 
about 1200km/s [600km/s] was covered by 512 channels with 
a width of 2.6km/s [1.6km/s] at 3mm [1mm]. The spacing be- 
tween individual grid points was 11", i.e. half the size of the 
3mm beam. The integration time per scan was usually 4min, 
and both polarizations were simultaneously observed for each 
frequency. Typical system temperatures during the observations 
were 350 K and 630 K for the 3mm and 1mm receivers. The 
data reduction was done using the GILDAS/CLASS software 
package. Each scan was inspected and those few with extremely 
high system temperatures or other instrumental effects were dis- 
carded. The baseline was corrected in the individual spectra by 
fitting a first order polynomial through channels outside the ex- 
pected line emission. After this correction all spectra for an in- 
dividual position were averaged together using a noise weight. 

2.3. Short spacing correction 

The 30 m observations were used to compute the short spac- 
ing correction (SSC) and recover the large-scale low-level flux. 



' http://www.iram.fr/lRAMFR/GILDAS 

- Robust weighting here means the weighting function W(u, v) in a 
iiv cell is set to a constant if the natural weight is larger than a given 
threshold, and W = 1 otherwise. 
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Table 1. Ovei-view of PdBI data 



Emission Line 


Weighting 


Beam Size 


PAO 


r.m.s. 






(" X ") 




(mjy/beam) 


''CO(l-O) 


natural 


2.56x1.70 


108 


2.8 




robust 


1.37x1.08 


114 


2.3 


'2CO(2-l) 


natural 


1.58x1.37 


87 


7.8 




robust 


0.64x0.50 


111 


8.0 



Notes: The beam sizes and PAs of the CLEAN beams for the PdBI data. 
The right most column gives the r.m.s. noise of the CLEANed PdBI- 
only data cubes. The data cubes are specified according to emission 
line and weighting. 



Table 2. Overview of PdBI+30m data 



Emission Line 


Weighting 


Beam Size 


PAO 


r.m.s. 






(" X ") 




(mJy/beam) 


''CO(l-O) 


natural 


3.11x2.59 


94 


2.5 




robust 


1.57x1.22 


112 


2.2 


'2CO(2-l) 


natural 


1.72x1.56 


81 


7.8 




robust 


0.69x0.55 


113 


7.8 



Notes: The beam sizes and PAs of the CLEAN beams for SSC (PdBI 
+ 30m) data. The right most column gives the r.m.s. noise of the 
CLEANed SSC (PdBI + 30m) data cubes. The data cubes are speci- 
fied according to emission line and weighting. 



The 30m observations were reprojected to the field center and 
frequency of the PdBI observations. The bandwidth coverage 
of the 30m observations was resampled to match the velocity 
axis of the interferometric observations. A combined data cube 
was produced using the task 'UV-short' in GILDAS. The sin- 
gle dish weight scaling factor was 6.55x10"^ for '^CO(l-O) and 
7.74x10-'* for i2cO(2-l). 

Two sets of final CLEANed data cubes were produced us- 
ing natural and robust weighting. The resolution of the natu- 
ral [robust] weighted '^CO(l-O) and '^€0(2-1) data is given in 
Table|2] The data cubes have 512x512 pixels, with a pixel scale 
of 0.25"/pixel [O.lO'Vpixel] and velocity bins of 10 [5] kms"' 
for the '2cO(l-0) ['^€0(2-1)] data. The r.m.s. noise per chan- 
nel in the SSC '^CO(l-O) observations is 2.5 mJy/beam and 
2.2 mJy/beam with natural and robust weighting, respectively. 
For the SSC '^CO(2-l) observations these values are both 7.8 
mJy/beam. CLEANing was done down to the 2cr noise level with 
the assistance of individual polygons defined for each channel 
with line emission present. 

The noise values of the SSC data cubes are, with the ex- 
ception of the '^CO(2-l) natural weighted data cube, below the 
noise levels of the PdBI only data. The beam sizes increase 
slightly (~ 8%) due to the added short spacings. 

2.4. HST/NICMOS observations 

We retrieved from the HST archive the NICMOS FllOW and 
F160W images of NGC 6951. The reduction was carried out us- 
ing the "best" calibration files, and the van der Marel algorithm 
(e.g. lBoker et al.ll 19991) to reduce the "pedestal" effect. Sky val- 
ues were assumed to be zero, which is ge nerally a go od assump- 
tion for these kinds of NICMOS images (iHunt & Ma lkan 2004). 
The images were calibrated, converted to magnitude scale, and 
subtracted to obtain a [Fl 10W]-[F160W]a; J-H color map. This 
color image and the starburst ring it reveals will be discussed in 
Sect. 3.1. 



3. Properties of the CO-emitting gas 

3.1. Morphology and H2 masses 



Table 3. Properties of NGC 6951 


Parameter 


Value 


Reference 


Type 


SAB(rs)bc 


(1) 


Nuclear Activity 


LINER/Seyfert 2 


(2) 


Dynamical Center 


(Locus Radio Continuum) 




RA (J2000) 


20*37"' 14. 123' 


(3,4) 


Dec (J2000) 


66°06'20.09" 


(3,4) 


Inclination Angle 


46° 


(5) 


Position Angle 


138° 


(5) 


Adopted Distance 


24.1 Mpc 


(6) 



References: (1) de Vaucouleurs et al.' (1991'), ( 2) iPerez et al.' ('2000^ 
(3) LE DA, KriDS et al. (2007), (4) Saikia et ah i2002h . (5) .Haan et S] 
J2009h . (6) ITullvllT988ft 



Table 4. CO line fluxes 



Component 


CO(l-O) 


CO(2-l) 




natural 


robust 


natural 


robust 




Sco 


Sco 


Sco 


Sco 


gas lane N 


23.3 


2.3 






gas lane S 


10.3 








spiral north 


142.4 


123.9 


280.3 


151.8 


spiral south 


109.9 


94.2 


202.8 


123.3 


central 


10.3 


6.1 


22.8 


0.9 


total 


296.2 


226.5 


505.9 


276.0 



Notes: Integrated line fluxes (Sco) for different components of the ob- 
served CO morphology. CO fluxes are given in ly km s ' . 



Table 5. M//, masses 



Component 


CO(l-O) 


CO(2-l) 




natural 


robust 


natural 


robust 










M„, 


gas lane N 


1.2 


0.11 






gas lane S 


0.5 








spiral north 


7.1 


6.2 


2.8 


1.5 


spiral south 


5.5 


4.7 


2.0 


1.2 


central 


0.5 


0.3 


0.2 


0.01 


total 


15.2 


11.6 


5.0 


2.7 



Notes: H2 masses (M/y,) for different components of the observed CO 
morphology. H2 masses in 10** Mq. For the C0(2-l) derived masses, 
we assume Ico(2-i)/Ico(i-0) = 0.8. 

The intensity maps of the CO(l-O) and CO(2-l) observations 
presented in Fig. [T]have been constructed f rom the CLEANed 
data cubes using the software 

GIPSM3 ( 

va n der Hulst et al.l 

1992). These zeroth moment maps are computed as the pixel- 
wise sum of emission above a fixed threshold. Here we chose 
the 3cr noise level. Spurious signals are filtered out by impos- 
ing the constraint that the emission be above this threshold in 
at least two consecutive velocity channels. All signals that sat- 
isfy these requirements are added to the intensity map. To the 

^ Groningen Image Processing SYstem 
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^ 30 1 



50 



40 



natural ^--g 

^^co(i-o.); 



spiral nortb 



robust 
^^CO(l-O) 



20 



S 10 



SB" 60' 



550 




4- ■■ 




30 



25 



g 66° 6 20" 



lar 16' 14" 12' 20''37'^10" 

Right Ascension (2000) 



natural 
^^C0(2-1) 




■11 



15 



10 



spiral south 



gas lanes 



I a' 16' 14" 12" 20^37"! 0' 

Right Ascension (200O) 



robust 

^^C0(2-1) 

' ''"SI 



16" 



15° 20''37"14' 
Riglit Ascension (2000) 



13" 



16" 



15" 20^37" 14' 

Right Ascension (2000} 



13* 



Fig. 1. Top: Integrated '2cO(l-0) emission in natural (a) and robust (b) weighting for the SSC data. The zeroth moment maps have 
been primary beam corrected. The '-CO(l-O) emission has been integrated from -220 km s ' to 250 km s"' . Contours run from Sa- 
in steps of lOcr. The r.m.s. value is: a) Icr - 0.025 Jy beam"' km s"' and b) Icr = 0.022 Jy beam"' km s"'. Bottom: Integrated 
'■^CO(2-l) emission in natural (c) and robust (d) weighting for the SSC data. The '^CO(2-l) emission has been integrated from -200 
km s"' to 200 km s"'. Contours run from c) 5cr in steps of 20o-, d) 3cr in steps of lOo". The r.m.s. value is: Icr - 0.039 Jy beam"' 
km s"' for both maps. The red cross indicates the position of the dynamical center in all panels (Table[3]l. The beam sizes are shown 
in the lower left comers and correspond to the values listed in Table |2] 



CO(l-O) emission zeroth moment maps we have applied a pri- 
mary beam correction. 

The intensity distribution of the observed CO(l-O) and 
CO(2-l) emitting gas has all the components we would expect 
from a gas distribution driven by a large-scale stellar bar. The 
morphology of the CO(l-O) emission map (Fig. \T\ top, left) 
from larger radii inward is comprised of the following compo- 



nents. At the edge of our field we see a spatially unresolved 
resonance ring at a radius of ~ 30" (=3.5 kpc). This is most 
likely at the lo cation of the ultra harmonic 4: 1 resonance inside 
corotation (e.g. lSakamoto et ani2000) . This ring is connected to 
straight gas lanes, with an approximately horizontal (east-west) 
orientation. Their position corresponds well to the orientation 
of the large-scale stellar bar, which has a position angle of 84° 
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25 



o 
o 
o 



a 

Q 



66" 6 20 




20"37"'14^ 
Right Ascension (2000) 

Fig. 2. '^CO(2-l) robust weighted intensity map (contours from 
3cr in steps of lOcr steps, Icr - 0.039 Jy beam"' km s"'); over- 
layed on a J-H color map (greyscale) from HST. White corre- 
sponds to higher extinction. 



dMulchaev et al.l 19971) . The gas lanes also coi ncide with the dust 
lanes at the leading edges of the stellar bar (iPerez et al.ll2000l 
Fig. 4a of their paper). In the robust weighted map of the CO(l- 
0) emission (Fig.[T] top, right), which is more sensitive to emis- 
sion on smaller scales, the straight gas lanes are no longer de- 
tectable. This indicates that the straight gas lanes consist of more 
diffuse, low intensity gas. 

Going further inward, we find the straight gas lanes connect 
to a 'twin p eaks' morpho l ogy an d a spiral pattern. This was also 
observed bv lKohno et al.l (Il999h . It is likely that the 'twin peaks' 
are due to the crowding of gas streamlines in a baiTed potential 
(jKennev et al. 1992). The orbit crowding leads to a buildup of 
gas at those locations where orbits change from xi to X2, i.e. 
with orientation along the bar major axis changing to orientation 
along the bar minor axis for elliptical orbits. The peaks inside 
the spiral arms have a distance of ~ 6" (- 700 pc) from the nu- 
cleus. Slightly within that radius ( ~ 5", 580 pc) a circumnuclear 
SB ri ng has been detected i n Ha ([Perez et al.ll2000l: iRozas et^ 
12002 ') and radio emission ( Saikiaet al.ll2002 '). 

The natural and robust weighted intensity CO(2-l) maps 
(Fig- El bottom) show a distribution very similar to the one seen 
in the CO(l-O) maps. Some of the differences, however, stem 
from the smaller field-of-view (FOV) of the CO(2-l) data. The 
outer ring at the edge of the CO(l-O) FOV is not visible, nor 
are the straight gas lanes. Due to the higher resolution, the two 
prominent peaks seen in the CO(l-O) maps now break up into 
multiple maxima. In the natural weighted map the northern peak 
is now joined by a second one slightly (0 ~30°) offset to the 
west. The southern spiral structure shows an elongated ridge 
with a far less distinct center. In the natural weighted CO(2-l) 
map the peaks and ridge are still unresolved structures. 

We start to resolve two tightly wound gas spiral arms, as 
well as peaks of dense gas within the ring only in our highest 
resolution robust weighted CO(2-l) maps. In the higher density 
tracer HCN(l-O) a similar effect is seen (Ki-ips et al. 2003). The 
peaks found in that high resolution image, although barely re- 
solved, sit almost 90° down the spiral arms. This seems to indi- 



CO 
>^ 

> 
> 




+ 10 +5 -5 
Rel. offset (arcsec) 



10 



Fig. 3. Position-velocity diagram along a PA of 1 13° through the 
dynamical center. Contours are at 2cr, in Icr steps, with the -2cr 
and - Icr contours given in grey. A gas bridge is visible between 
the northern part of the ring (lower left here) and the nucleus 
(denoted with a cross). 



30 



o 
o 
o 

g 66° 620" 








16^ 



20"37"'ir 
Right Ascension (2000) 



Fig. 4. Emission line ratio CO(2-1)/CO(1-0) in temperature 
units. This ratio was computed using the natural weighted CO(2- 
1) and the robust weighted CO(l-O) data cubes, smoothed to the 
resolution of the former. The high ratio south-east of the nucleus 
inside the ring is not significant. 
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cate some density gradient in the gas ring, starting at the north- 
ern and southern peaks and going downstream. Fig. |2]shows the 
CO(2-l) integrated intensity map as contours superimposed on 
the NICMOS J -H color map (Sect.|Z4li. The figure shows that 
the CO spiral arms coincides nicely with the SB ring seen in the 
red J - H color. The extinction. Ay, can be derived from this 
color image assuming aJ-H color of ~0.7 for a normal (unred- 
dene d) composite stellar p opulation, and (Aj - AeVAy = 0.092 
(as in lCardelh et"!!] d 19891) . with TJy = 3.1). The Ay's so derived 
for the SB ring are quite red, ranging from ^1 to more than 6 in 
the dustiest, most compact regions. 

Finally, at the dynamical center (Table |3]l we find barely re- 
solved molecular g as emission (Mx 10^ Mp^) a t the 3cr level, as 
already reported by iGarcia-Burillo et al Central molec- 

ular gas emission has also been detected in HCN at a much more 
prominent level (Krips et al. 200 71). The g as bridge between the 
ring and the nucleus, identified by Garcia-Burillo et aD (l2005h as 
lying at P. A. 134°, is also observed in our data cubes. Here we 
find the bridge to be present between PAs of 78° and 158°, but 
the most intense emission is along a PA 113°, as shown in the 
position- velocity diagram in Fig. |3] 

The morphology discussed here was for the most part also 
observed in the previously published PdBI-only maps. The ex- 
ception is the spatially unresolved resonance ring at the edge 
of our field, which was not seen at all. The straight gas lanes 
are more prominent in the SSC maps presented here. The conse- 
quences for the gravitational torques will be discussed in Section 
lO 

In Tables |4] and |5] we list the integrated CO line fluxes and 
the derived H2 masses for the different components described 
above. The regions over which the flux has been measured are 
indicated in Fig.[TJ). These values have been derived with a CO- 
to-H2 convers ion factor Xco of 2.2 x 10^*' cm"^ [K km s"']"' 
(ISolomon & B arrett 1991) and have not been coiTected for he- 
lium abundance. As the CO line ratio (Sect. 13.2b is fairly con- 
stant in the center, the assumption of a single conversion factor 
seems valid for our purpose. There are claims that the CO-to- 
H2 conversion factor is a factor 3-4 lower in gas-rich centres of 
spiral galaxies (e.g. WeiB et al. 2001). A lower conversion factor 
would lower our mass estimates coiTespondingly. However, for 
consistency with other papers in this series, we use the Galactic 
value. The total integrated '^CO(l-O) flux within a 40" field of 
view is 314 + 55 Jy km s The '^00(1-0) flux we find, is 
si milar to the measure ments of 334 + 12 Jy km s"' within 65" 
bv lKohno et alj (i 19991 obtained with the Nobeyama Millimeter 
Ar ray and the 45m tele scope) and 350 + 41 Jy km s ' within 45" 
bv lYoungetaLl (119951 with the FCRAO single dish telescope). 
The corresponding H2 mass is 1.6 x lO** Mq. For the '2CO(2-l) 
emission we measure an integrated flux of 452 + 65 Jy km s"' 
within a field of view of 14", which corresponds to an Mh, of 
4.6 X 10** Mq if we assume Ico(2-i)/Ico(i-o) - 0.8 (see Sect. 13.2b . 
As can be seen from Table |5] the total measured mass in CO is 
completely contained in the components mentioned. The areas 
within which flux is measured are the same for natural and ro- 
bust weighted maps. In both CO(l-O) and CO(2-l) lines we see 
that the natural weighted maps are more sensitive to low-level 
diffiise components of the CO gas. The areas over which the flux 
of the components has been measured, are marked in Fig. [T] 

3.2. Line ratio 

The CO(2-1)/CO(1-0) line ratio has been used to convert the ob- 
served CO(2-l) flux into equivalent CO masses in Section ITTI 
We have derived the line ratio (Fig. |4) in the following manner. 



Both maps, CO(l-O) and CO(2-l), have been short-spacing cor- 
rected and sample down to the same minimum uv radius. The 
robust weighted CO(l-O) data cube was smoothed to the reso- 
lution of the CO(2-l) natural weighted cube. Then the zeroth 
moment map of this smoothed CO(l-O) cube was constructed 
as before. The CO(2-l) zeroth moment map was regiidded to 
the pixel scale of the CO(l-O) map (from 0.10" to 0.25"). The 
fluxes of both maps were converted into temperature and the ra- 
tio taken. 

The ratio is almost constant along the spiral arms, with the 
map displaying an average ratio of 0.8. The high ratios found at 
the edges of the ring are insignificant due to low S/N in the CO 
maps. 

3.3. Kinematics 

Evidence of the dominant influence of the large-scale stellar bar 
is visible in the velocity field. CO(l-O) line emission has been 
detected in the velocity range of -220 km s"' to 250 km s"' 
relative to the heliocentric velocity of 1425 km s ' in the natu- 
ral weighted data cube. For the CO(2-l) emission, the velocity 
range is slightly smaller; from -200 km s"' to 195 km s"'. That 
is in part due to the higher rms noise in this data cube, which 
affects detection of the signal in the channel maps at the highest 
relative velocity offsets. All channel maps with significant emis- 
sion are shown in Figs. |5]and|6] As the velocity increases, the 
line emission shifts from the north-west to the south-east. This 
is consistent with the major kinematic axis of this galaxy, mea- 
sured by Haan et al. (2009) as 138°using HI data (TableO. The 
channels close to systemic velocity (-50 to 100 km s"') show 
two maxima, from both sides of the CO ring, as well as extended 
arms from the bar-driven straight gas lanes. 

The CO( 1 -0) velocity map is shown in the top left panel of 
Fig-El The iso-velocity contours in the center are almost perpen- 
dicular to the major kinematic axis. As the radius increases the 
iso-velocity contours bend, forming the 'S'-shape distinctive of 
velocity fields dominated by large-scale bars. 

The dispersion (see the velocity dispersion map in the right 
panel of Fig.|7]i reaches values of up to 70 km s"^ These values 
are high for a gas disk and have led us to further investigate the 
kinematics of the CO. The data cubes show that the observed 
CO line emission arises from two distinct components in veloc- 
ity space at several positions, most prominently where the dis- 
persion maps show high values. This complex velocity structure 
means that we are not seeing a truly high local velocity disper- 
sion in a single component, but rather the projection of two ve- 
locity components within the same beam. The assumption of a 
single component when we determine the velocity dispersion is 
clearly wrong for some positions. We see these double peaks in 
both CO(l-O) cubes and the natural weighted CO(2-l) cube. 

In order to quantify the double emission peaks in detail, we 
fitted double Gaussians along the spectral axis at each spatial 
pixel of the CO(2-l) natural weighted data cube (Fig. [8]l. The 
choice of this cube was made based on the higher spectral res- 
olution of the (2-1) cube with respect to the CO(l-O) cubes. We 
used the function 'XGAUFIT' from GIPSY for the fitting. At 
the spatial pixels where we have a double peak, we separate the 
two components based on the large-scale bar model we derive 
in Sect. 14.21 In Fig. ^{middle, left) we show the central value of 
the fitted Gaussian for pixels requiring only a single Gaussian fit, 
and the central value of the Gaussian fit closest to the bar model 
for the pixels with a double Gaussian fit (component 'Vl'). This 
results in a good representation of the velocity field of the disk 
of the galaxy. In Fig.^(bottom, left) we present the central val- 
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ues of the second Gaussian (component 'V2'). For the most part, 
the double peaks are present in the northern spiral arm/ring re- 
gion, with the exception of the complex in the south-south-west 
close to the nucleus. The velocity difference, Av, between the 
two components is between 40 km s"' and 120 km s"' . 

Most of the two components might be interpreted in terms 
of xi and X2 orbits, xi orbits are parallel to the large-scale bar 
(building up the straight gas lanes), and X2 orbits are perpendic- 
ular to the bar (supporting the ring). At the resolution of the data 
cube (beam size 1.72" by 1.56"), the two orbit families might 
blend together. We find double emission peaks in the CO(l-O) 
cubes and the natural weighted CO(2-l) cube, where we are un- 
able to resolve the end of the one spiral arm from the other spiral 
arm. Double velocity components would be the natural conse- 
quence of the emission of CO gas on the two orbits families 
being blended together spatially in our data cubes. In compar- 
ison, in the robust weighted CO(2-l) cube we do not find ev- 
idence of double peaks in velocity. We find little difference in 
the width (velocity dispersion) of the Gaussians (Fig. Q mid- 
dle/bottom, right), except in the northern region. There we see 
a narrower component, connected to 'Vl', with a dispersion of 
about 30 km s"', and a wider component, connected to 'V2', 
with a value of around 40-50 km s '. This region is where we 
have the northern peak of the 'twin peaks' morphology, which 
persists even in the high resolution CO(2-l) intensity map (in 
the south we see a ridge). The high dispersion here might be ex- 
plained by the gas being shocked, which leads to more real, local 
turbulence/dispersion. 

The resolution argument does not explain the double emis- 
sion peaks found inside the spiral arms, south-south-west of the 
nucleus. In this region, based on the CO(2-l) map, we expect 
very little emission, let alone emission arising from two distinct 
(Av ~ 80 km s"') components. The velocity of the first compo- 
nent agrees well with what is expected from the large-scale bar 
velocity field. The second component's central velocity (~ 70-90 
km s"' ) seems to connect it with the southeastern part of the ring. 
There is little difference in dispersion between the two compo- 
nents. The strength of the second Gaussian is nearly equal to the 
first component. This double emission region is not spatially co- 
incident with the gas bridge mentioned earlier. In Section l672l we 
will discuss the significance of both the bridge and this second 
kinematic component in the nuclear region for nuclear fueling. 




Fig. 8. Six example spectra from the robustly weighted CO(2-l) 
data cube where we fitted a double Gaussian to the velocity axis. 
The location of each spectrum is indicated with red lines in the 
intensity map (middle; same area as Fig.[T]bottom,left panel). 



4. Kinematic modeling 

4.1. Motivation 

NGC6951 has a prominent large-sca le stellar bar with a ra dius 
of 3.0 kpc at a position angle of 84° ( Mulchaev et al.ll 19971) . We 
suggested in Sect. |3]that the gas morphology and kinematics in 
NGC 695 1 can be well explained in terms of this barred poten- 
tial. The spiral pattern and the high intensity twin peak morphol- 
ogy are both consistent with gas streamlines under the influence 
of a large-scale bar. 

To further investigate this claim, we model our CO ob- 
servations with a gravitational potential where the only non- 
axisymmetric component is a large-scale bar. This cannot be 
done from observations directly since this would assume ex- 
act prior knowledge of the non-axisymmetric components of the 
gravitational potential. At the same time this modeling will also 
allow us to derive gravitational torques in this region for an in- 
dependent comparison wit h th e gravitational tor que results of 
iGarcfa-Burillo et al.l (l2005l) and lHaan et aLl(l2009h (Sect.Ol- 

4.2. Description of the model 

We modeled the inner 5 kpc (diameter) of the observed CO(l-O) 
gas distribution and its underlying gravitational potential using 
the DALIA modeling software dBoone et a l. 2006) that was al- 
ready used to model the NUGA target NGC 4569 bv .Boone et al.l 
(I2007h . 

The tool constructs a kinematic model of gas particles in the 
presence of a barred potential. The potential is built up from two 
components, an axisymmetric component representing the disk, 
and a weak m=2 perturbation: the 'bar' . The axisymmetric com- 
ponent has a logarithmic shape and is defined by a characteristic 
length (r^,) and velocity (Vp). 



(1) 



The bar has the same radial profile as the axisymmetric com- 
ponent, but is tapered by a sine in azimuth and can be oriented 
in any preferred direction relative to the galaxy's kinematic ma- 
jor axis (Table |6] perturbation azimuth). It is further defined by 
the relative amplitude of the perturbation with respect to the ax- 
isymmetric component (i.e. the bar strength, e), and its pattern 
speed. 



^b{r, (/>) = e<l)()cos(20) 



(2) 



The total kinematic model potential can therefore be described 
as follows: 



(1 + ecos(20))<l)o 



(3) 



The model is populated with a distribution of gas particles, 
described in terms of the radially varying parameters column 
density, velocity dispersion and scale height. Two dissipation 
terms, one acting radially and one acting azimuthally, are set to 
reproduce the dissipative behavior of the gas particles. Finally 
the model can be inclined and positioned corresponding to the 
observations, and a rotation sense is set. 

This model h as several limitations as dis cussed by 
Boone et al 1(l200l. Firs t, it is singular- at corotation. iButa et all 
(l2003h and lBlocketal.1 (|2004 place NGC6951's bar radius at 
3.0kpc and indeed, as can be seen from Fig. |9](to/7, right), our 
model does not extend to corotation. Second, closed orbits are 
computed, which are not self-consistent with the inclusion of 



8 T.P.R. van der Laan et al.: Molecular gas in NUclei of GAlaxies (NUGA) 




_Ll I I L 



Fig. 5. Channel maps of the naturally weighted '^CO(l-O) data cube. The size of each channel map is 44" by 44". The contours 
are at -3cr, -2cr, 3cr, Scr, lOcr, IScr and 25cr, with Icr = 2.5 mJy/beam. The velocity relative to systemic velocity of the galaxy 
(Vsys = 1425 km s ') is indicated in the upper left corner. The dynamical center is indicated by a cross in each channel map. The 
synthesized beam is given in the lower left comer of each channel and the dirty beam is shown in the lower right panel. 



two dissipative terms. Therefore, the locations of resonances ob- 
tained from this model should be taken as a first-order estimate 
only. 

All parameters in the model can in principle be freely cho- 
sen. Our interest here is in obtaining a good set of parameters 
that reproduces well the observed CO gas distribution. For this 
reason, we adopted first guesses for the model parameters based 
on known values of the system and subsequently slightly varied 
them so that the model better fits our observations. The radial gas 
mass distribution (Table|7]i was deduced from the CO(l-O) emis- 
sion maps (Fig.[Tl top) and can be scaled with the total observed 
H2 mass. For the velocity dispersion we take a single value, 30 
km s"', which is reasonable for these radii and th e assumption 
of a pressure supported gas disk (iJogee et al.ll2005b . The rotation 
sense is set to clockwise, as the dust lanes along the bar are as- 



sumed to be leading. The dissipation rates were initially left at 
default values and only changed after the other parameters had 
been fixed. Table |6] lists the values of the potential parameters 
with which the model best reproduces our observations. The best 
fit was estimated by comparing the modelled channel maps (Fig. 
[Tot and intensity map [Fig.^bottom, left) against their observed 
counterparts by eye. In Appendix A we show the agreement of 
our best fit model with the CO observations at several different 
positions in the data cubes. 

4.3. Best fit model 

Figure |9] shows the orbits {top, left) and rotation curve {top, 
right) for our best model. The orbit pattern shows the change 
from xi to X2 orbits, the main periodic orbits in a bar poten- 
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Fig. 6. Channel maps of the naturally weighted '^CO(2-l) data cube. The size of each channel map is 25" by 25". The contours are 
at -3cr, -2cr, 3cr, 5cr, lOcr, 15cr and 25cr, with Icr - 7.8 mJy/beam. All other notation is as in Fig.|5] 



tial, starting at r=1.5kpc. The existence of two ILRs, denoted 
ilLR (inner) and oILR (outer) in our model, can be seen from 
the double intersection of the Q - f curve with the pattern speed. 
Corotation is beyond 2.5 kpc, which means the model does not 
diverge in the range we are studying here. 

The presence of two ILRs in NGC 695 1 has a lso been found 
by both Rozas et aL (.20021) and lPerez et al.l (l2000l ). who put them 
at 180pc and 900 pc, and ISOpc and llOOpc, respectively. Our 
CO observations cannot be fit by the model without 2 ILRs. 
Their locations are determined by the free parameters of the 
model; the pattern speed and the shape of the potential (defined 
by its characteristic length and velocity). From the resonance di- 
agram in Fig.|9]we find our ILRs to be at 160pc and 1150pc. 
These radii are close to the estimates reported by the earlier pa- 
pers. The orientation of the bar perturbation (13 0°) in the model 
is also close to the value for the bar reported bv lMulchaev et al] 
(Il997h . This degree of consistency shows that the model is a 



good representation of the observations, and thus we conclude 
that the CO gas kinematics are indeed dominated by the large- 
scale stellar bar 

The final model moment maps and channel maps are pre- 
sented in Figs. |9] (bottom) and (TO] The model has been com- 
puted at the resolution of the CO(l-O) data. A comparison of 
Fig. [To] with the observations (Fig. |5]l shows good agreement. 
The model reproduces the straight gas lanes as well as the two 
peaks of emission in the ring in the individual channels from -30 
km s ' to 40 km s '. The channels with emission in the model 
are only slightly fewer in number than in the observed channel 
maps. In the integrated emission map of Fig.^(bottom, left), the 
two spiral arms and the twin peaks are also well reproduced in 
this model. It is very encouraging that this bar potential model 
with its simple kinematics and intrinsic limitations fits the ob- 
servations so well. Therefore, we are confident when we use the 
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Fig. 7. Top left: First moment map of the CO(l-O) emission. Contours: -200 to 175 km s"' with 50 km s"' steps, km s"' thick 
gray contour. Kinematic major axis (white dashed line at PA = 138°), orientation of the large-scale bar (dashed-dotted line, PA 
= 84°). Top right: Second moment map of the '^CO(l-O) emission. Contours: to 70 km s"' in 10 km s"' steps. The dynamical 
center is marked by a black cross. Middle/Bottom panels: Decomposed velocity field of the CO(2-l) natural weighted line emission. 
Black spots correspond to blanks due to bad fits at some spatial pixels. Middle left: The central value of the Gaussian fitted for each 
spatial pixel. If two Gaussians where fitted, the value of the Gaussian closest to the large-scale bar model velocity is given. The 
white contours represent the CO(2-l) intensity map (5cr in 40cr steps). Bottom left: For the spatial pixels where a double Gaussian 
was fitted, we plot the central value of that second Gaussian here, i.e. the value further away from the bar model velocity. The size 
and PA of the nuclear stellar oval and PA of the gasbiidge, that will be discussed in Sect. 16.21 are highlighted with a white ellipse 
and dashed lines. Middle/Bottom right: Velocity dispersion corresponding to respective velocity components. The dispersion is the 
Gaussian fitted cr at each position. Contours from to 70 km s The region with significant differences in velocity dispersion 
(discussed in Sect. l3.3l l is highlighted with a white ellipse. 
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model velocity field to separate the two velocity components we 
detected in the data cubes. 

In the channels, -30 km s ' to 40 km s the orientation 
of the two emission peaks is slightly more north-south in the 
model channel maps than in the observed one. From Fig. Qwe 
already saw that when we decompose the velocity field with two 
Gaussian fits, the velocity contours of the 'bar' -component also 
become more north-south oriented. 

In Sect. 13.11 we reported the location of higher density gas 
traced by CO(2-l) and HCN further downstream in the spiral 
arms than what is seen in the CO(l-O) line, which is mainly trac- 
ing lower-density gas. High emission peaks in both the CO(2-l) 
and HCN gas distribution occur away from the 'twin peaks' that 
dominate the CO(l-O) distribution. We attemped to reproduce 
this pattern by also computing a best model at the resolution of 
the CO(2-l) data. The model gas distribution still agrees with the 
observations in its main characteristics; spiral arms and the twin 
peaks. However, the details, such as the southern intensity ridge 
and the second peak in the northern spiral, are no longer well 
reproduced. We point out that this model therefore reproduces 
a morphology closer to the CO(l-O) emission than either of the 
higher gas density maps mentioned above. A plausible explana- 
tion could be that the CO(2-l) and HCN distributions reflect the 
evolution of the molecular gas inside the circumnuclear ring it- 
self, after the large-scale bar has brought the molecular gas there. 
As the molecular gas streams within the ring, cloud collapse can 
occur and/or continue, leading to higher molecular gas densities 
away from the contact points between ring and gas lanes. 

Table 6. Bar potential model parameters 




16° 15° 20^37'"! 4' 13" 12" 
Right Ascension (2000) 



16" 15"20''37'^14"13'' 12" 
Right Ascension (2000) 



Parameter 


Value 




Initial Guess from data 


Inclination (deg) 


46 


46 


Position Angle asc. node (deg) 


318 


318 


Char. Length (pc) [r,,] 


230 




Char. Velocity (km s"') [Vp] 


200 




Perturbation Azimuth (deg) 


130 


84 


Sense of Rotation 


Clockwise 


Clockwise 


Bar Strength [e] 


-0.035 




Pattern Speed (km s"' kpc"') 


50 




Rad. Diss. Rate (km s"' pc"') 


0.000 


0.005 


Az. Diss. Rate (km s"' pc~') 


0.020 


0.002 



Notes: The effect of errors in these parameters on our results is dis- 
cussed at the end of Sect. [5] 



5. Model derived results 

The model described in Sect. 14.21 can also be used as the basis 
for a gravitational torque computation. The model's large-scale 
bar potential is known at each position, since the model is purely 
analytical. This makes it easy to derive the gravitational torque 
T(r, 0) on the gas particl es at each point, following the equation 
from lBoone et all boOTi) : 



f(r,4>)^ -ev],ln{\ + —)sin{24>) (4) 
'"/' 

Here e is the bar strength, Vp the characteristic velocity and 
the characteristic length as set in the model. ^ is the azimuthal 
angle with respect to the bar position angle. 

Since the orbits for the gas particles are known, we define the 
net torque at a given radius as the time-average integrated torque 



Fig. 9. Top, left: Deprojected orbit pattern in the inner 3.0 kpc of 
the bar model. Top, right: Resonance diagram. The rotation fre- 
quency Q (red), - f (blue), + f (cyan) and the bar pattern 
speed (yellow dashed). The ilLR at 160 pc, oILR at 1 150pc and 
CR at 4000 pc (outside the displayed range) are indicated with 
black dots. Overlayed is the rotation curve, which approaches 
200 km s"'. The vertical dashed line indicates the outer edge 
of the models radial mass distribution. Bottom, left: Integrated 
emission of the bar model with contours of CO(l-O) intensity 
map overlaid. The model emission has been computed at the 
CO(l-O) resolution. Contours run from 5cr in steps of 20o- (Icr 
— 0.025 Jy beam"' km s"'). The red cross indicates the position 
of the center. Bottom, right: Velocity field of the bar potential 
model with contours of the CO(l-O) velocity field. Contours run 
from -200 to 200 km s"' in 25 km s"' steps. 



over the orbit with this characteristic radius. We recall that the 
characteristic radius of an orbit is the radius this orbit would have 
without a bar (the orbit would then be circular). As mentioned in 
Sect. 14.21 this approach is not completely self-consistent because 
when there is angular momentum change the orbits cannot be 
perfectly closed. Therefore, the gravitational torques computed 
here should be considered as a first-order approximation to the 
case in which the orbits are nearly closed. This is a valid approx- 
imation if the angular momentum loss is small. Using the simple 
argument that ^ is of order where v is the characteristic ve- 
locity of the potential v=200 km s"', we can see that the relative 
amplitude of the angular momentum loss per orbit is of order 
~ 0.0125 - 0.05. Therefore, it takes at least 20 orbital times for 
the gas to fully lose its angular momentum. 

Figure [TT](/e/0 shows the net gravity torque for each orbit 
in the model. On the horizontal axis we plot the characteristic 
radius of each orbit. The net gravity torque is negative over the 
entire radial range studied here. From 1.9 kpc down to 1.0 kpc 
the torque stays somewhat constant at values of -2000(km/s)^. 
There is a slight upturn at larger radii. This is the region of 
the straight gas lanes. The large value indicates that the bar is 
easily able to transport gas inward to the gas spiral arms. The 
torques then increase at an almost constant rate down to a ra- 
dius of 400 pc. We can understand this from Fig. ^{top, left): 



12 T.P.R. van der Laan et al: Molecular gas in NUclei of GAlaxies (NUGA) 

Table 7. Bar potential model radial density parameters 



Radius (pc) 


20 


450 


500 


600 


700 


800 


900 


1200 


1500 


2500 


Col. Dens. (rel. scaling) 


0.25 


0.2 


0.6 


0.7 


0.6 


0.6 


0.3 


0.1 


0.1 


0.02 


Velocity Disp. (km s"') 


30 


30 


30 


30 


30 


30 


30 


30 


30 


30 


Scale Height (pc) 


10 


10 


10 


10 


10 


10 


10 


10 


10 


10 



the orbit position angle with respect to the bar stays almost con- 
stant (at its maximum value) and the orbit elongation decreases 
with decreasing radius. A second, smaller, minimum is present 
at r=150pc, at the ilLR. The net gravity torque becomes zero 
within the 50 pc radius. A close inspection at the model orbits in 
the inner 300 pc shows that the orbit orientation changes again 
from JIC2 to jci. 

The mass accumulation rate that we can calculate from the 
gravitational torques also depends, in contrast to the gravita- 
tional torque per se, no longer solely on the analytical poten- 
tial, but also on the radial gas mass distribution in the model. 
Fig.[TT](r/g/if) shows the mass accumulation rate for this model. 
The variation in the mass accumulation rate indicates that gas 
must accumulate on some orbits and be depleted from others. 
That the gravitational torques drive the gas inward is also evi- 
dent from this curve. At radii larger than 800 pc gas is being de- 
pleted and the accumulation rate is negative. These are the radii 
where the gas is located in the straight gas lanes. From 400 pc to 
800 pc radius, the accumulation rate becomes positive. These are 
the radii at which we find the gas spkal arms and the gas ring, 
which are separated into two peaks in the figure. In the radius 
range 400-1 000 pc we already saw from the torque plot that the 
torques decrease with decreasing radius. As the radius decreases 
the torques are less strong, so less gas is being moved inwards at 
each radial bin and the gas can accumulate. Inside 400 pc we see 
a sharp positive peak, at 150pc. This is the radius at which we 
found a local minimum in the gravitional torques. The feature 
is confined to only a very limited radial range, but indicates that 
some accumulation may still happening inside the ring. 

Mass accumulation within 1.5 kpc in radius, i.e. integrated 
over this radial range, is occuring at a rate of ~+2.0 M0yr"' 
and mostly in the two spiral arms. From our observations we 
know that the total amount of molecular gas within this radius is 
about 1.6xlO^M0. Therefore, at our derived accumulation rate, 
all matter in the spiral arms must have been brought there within 
the last 0.8 Gyr However this is a lower limit to the age of the 
ring, since it assumes that no gas has been turned into stars and 
that there is always enough gas available at larger radii to be 
driven inward. 

Errors in the gravitational torque computation can come 
from six model parameters; the characteristic length and veloc- 
ity of the potential, the bar pattern speed, the bar strength and 
the dissipation factors. The first three influence the locations of 
the resonances as well as the values of the net gravity torques. 
A 10% difference in either of these parameters can cause up to 
a ~20% difference in the locations of the gravitational torque 
minima. The change in value is 5% for the characteristic length 
and 20% for the characteristic velocity and pattern speed. The 
net gravity torque scales linearly with bar strength. The dissi- 
pation factors only influence the strength of the gravitational 
torques. A 10% difference in these factors leads to a decrease of 
X! 25% in the net gravity torques. However, such changes to the 
model parameters lead to significantly worse models. Therefore, 
the uncertainties in the results we derive here for the net grav- 
ity torque and the mass accumulation rate are within the relative 
errors given here. 



6. Discussion 

6.1. Gravitational torques 

We can now compare different gravitational torque measure- 
ments w e obtain using either the ' torque map' method intro- 
duced by Garcfa-Burillo et ^0(120051 ; hereafter: GB05) (see also 
Ha an et al. 1 120091 ; hereafter: H09), and the analytical co mputa- 
t ion po ssible for the large-scale bar model following Boon e et"al] 
(l2007h . 

First, we investigate how the inclusion of the 30m data 
changes the results obtained with the 'torque map' method. 
Using the code PyPot of H09 we have derived the gravitational 
torques for our PdBI+30m data cube. PyPot uses NIR high- 
resolution images to evaluate the stellar potential. The gravita- 
tional torque curves for PdBI-only (dashed line) and PdBlH-30m 
(solid line) data are shown in Fig. [12] As can be seen from 
the figure, the inclusion of the 30m data leads to a change in 
the gravitational torque budget, especially in the region inside 
500 pc. The high positive torques derived at 400 pc in the PdBI- 
only maps disappear in the PdBI+30m result. This can be ex- 
plained by the larger amount of diffuse molecular gas that is re- 
covered in the PdBI+30m data, which lowers the high torque 
baiiier at -200 pc found by GB05 (see their Fig. 12) and H09. 
The spiral arms become more prominent with the inclusion of 
the 30m data, so they contribute more to the torque budget at 
larger radii, shifting the minimum to larger radii. This test shows 
that sampling the line emission on all spatial scales is important 
when computing gravitational torques from observations as the 
absolute values can vary significantly while the overall shape is 
preserved. 

Now, we turn to the comparison of the 'torque map' result 
with the analytical computation. We are interested in the torques 
as we want to compute the angular momentum change in the gas 
(whether the gas is in- or outflowing). In the analytical case, we 
can simply integrate the net torque over gas orbits from the bar 
model potential (which is a simplification of the true underlying 
potential). This integration gives the angular momentum loss of 
a gas particle in one orbit. In the case of 'torque maps' it is often 
argued that it is impossible to know the exact orbits from ob- 
servations. So, a statistical argument is used in the 'torque map' 
method. The assumption is that the observed gas distribution can 
be linked to the time spent by the gas along the orbits. The re- 
sulting torque is the gas column density weighted average of the 
local torques averaged per radial bin. 

There seem to be significant differences between the re- 
sults form both methods. In the 'torque map' result we find the 
strongest (negative) torques at a radius of ~650pc. For the ana- 
lytical computation the most negative torques are approximately 
at a radius of r~1200pc. As we saw from the intensity maps (Fig. 
[TJ, we observe very little gas in the straight gas lanes (r> 1.0 kpc). 
We know, however, that the large -scale stellar bar has a radius 
of 3.0 kpc dMulchaev et alj[T997l) . At 1 kpc the large-scale bai" 
still exerts significant torques on the gas, but the 'torque map' 
method being very sensitive to the actual observed gas distri- 
bution would not show that since there is so little gas, hence 
the earlier upturn. There is also a difference in the absolute val- 
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Fig. 10. Channel maps of our best model fit at the resolution of the CO( 1 -0) data. The FOV of each channel map is 44" by 44" . The 
velocity indicated in the upper left comer is relative to the central channel (0 km s ' ). 



ues of the strongest negative torques between the 'torque map' 
method and the analytical computation. This inequality is due 
to the difference between computing the net torque (model) and 
averaging the torques (maps). 

In the inner 300 pc, another difference is seen between the 
'torque map' method and the analytic computation. In our an- 
alytical computation we see a second small dip in the torques, 
while in the 'torque map' result we find positive torques. By con- 
struction, it is impossible in our model for the torques to change 
sign. Further, GB05 have found that there is a secondary grav- 
itational component, a nuclear stellar oval, that plays a role in 
the torque budget here. No such secondary component has been 
included in our model. 

Despite these differences, the comparison between our mass 
accumulation rate in Fig. [TT] and the mass accumulation rate 



found by H09 (see their Fig. 12) shows little difference. We find 
values not more than a factor 2 larger than theirs. This agree- 
ment shows that our large-scale bar potential model and the sub- 
sequent analytical computation give a good representation of the 
true gas transportation in this region, at least from the outer disk 
down to the circumnuclearring. 



6.2. CO gas flow inside ttie ring 

Another aim of this paper is to see if a large-scale bar can ex- 
plain the gas distribution and kinematics over the entire 1.5 kpc 
radial range. For all radii down to the ring the model repre- 
sents our obsei-vations very well. Inside the ring we found neg- 
ative torques in our model, but positive torques based on H09's 
'torque map' method. Does the influence of the large-scale stel- 
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Fig. 11. Left: Net gravity torque derived from the best fit model. Right: Gas accumulation rate as a function of radius. The orbital 
characteristic radius corresponds to the radius of circular orbits that a given gas particle would maintain if unperturbed by the 
model's large-scale bar. 
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Fig. 12. Gravitational torque curve deriv ed from th e PdB I- 
only data {{dashed line) reproduced from iHaan et all (l2009l )). 
Gravitational torque curve derived from the PdBI+30m data 
{solid line) using the same method. The addition of the 30m 
data has a large influence on the torque values, especially inside 
500 pc. 



lar bar extend further in? As a first test we looked for evidence 
of gas spi ral arms inside the ring, as predicted by the model of 
lEnglmaier & Shlosm an (2000) and Maciejewski (2004). In Fig. 
12] we see some indication of spiral structure in the dust. In the 
observed CO morphology we do not, nor do we fi nd such evi- 
dence when we decompo se the velocity field ( Fathi et aDl2005t 
Ivan de Ven & Fath ?'201CA. The latter method is sensitive to even 
a small arm/inter-arm density contrast. Based on this outcome, 
it seems that the CO gas is mainly on stable orbits in this region, 
which implies that it is not inflowing. However, we did find other 
evidence of CO gas that is potentially not on stable X2 orbits; the 
gas-bridge and the double CO emission peaks in the line profiles 
(Figs. SandEl Sect. Elll. 

The primary component of the double peak emission inside 
the ring (Fig.|7] middle right) shows velocities indicative of sta- 
ble orbits. The line-of-sight velocity of the 'second peak' CO 
component (located within the white ellipse in the bottom right 
panel of Fig. [T) is 70-90 km s"', which corresponds to a velocity 



difference from the primary CO component at these coordinates 
of Av ~ 80 km s"' . It is unlikely that we have a simple superposi- 
tion of two gas clouds, both on stable orbits. After separating the 
flux in this area into the two velocity components, we derive an 
H2 mass for this 'second peak' CO component of I.SxIO^Mq. 

To understand if this CO component might still be related 
to gas spiral arms induced by the large-scale bar, the gas flow 
direction of the 'second peak' component and whether it is lo- 
cated in the disk of the galaxy must be ascertained. Th e first issue 
is straightforward. As Storchi-Bergmann et al.l(l2007l have con- 
cluded the near side of the galaxy is to the southwest and the 
far side to the northeast. Thus, the positive velocity of this CO 
component, which is south-south-west of the nucleus, indicates 
streaming toward the nucleus if the component is in or in front of 
the disk, and outflow if it is behind the disk. Outflows are usually 
biconical. If this component is part of an outflow, seen behind the 
disk, we would expect to see a second blue-shifted component 
from the outflow in front of the disk, which is not the case. Nor 
do we really expect an outflow, since this galaxy does not have a 
strong AGN or starburst near its nucleus. It is thus unlikely that 
this CO component is located behind the disk and outflowing. 

If the component is inflowing, it is very unlikely that it has 
any significant scale height above the disk. NGC 695 1 is an iso- 
lated galaxy that has been dynamically undisturbed for at least 1 
Gyr. This history means that the origin of the CO component has 
to be internal. It could be that some intense star formation in the 
ring expelled CO from the ring. We do see a gap in the CO ring 
close to the component. However, the low velocity dispersion of 
the CO gas in the component makes it unlikely that the CO com- 
ponent was heated due to star formation, became unbound and 
was pushed high above the disk and is now falling back. Nor do 
we see evidence of a recent starforming event at that location. 

Thus, it seems most likely that the CO component is inflow- 
ing and located in the disk, bringing us back to the question of 
whether the component's motion can be explained by large-scale 
bar induced gas spiral arms, even though our initial tests were 
negative. Storchi-Bergmann et al. (2007) claim to see evidence 
of two such spiral arms in Ho' line emission. The location of 
their Ho- spiral arms does not match the location of our CO com- 
ponent however. Their derived velocity of + 20 km s ' is also 
at least a factor 3 smaller than that of our CO component. If 
there are really gas spiral arms in NGC 6951, this discrepancy is 
very unexpected. CO should be a better tracer of such gas spi- 
ral arms, since it is the dominant gas tracer in this region. So it 
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seems that the large-scale bar does not induce further gas inflow 
in this galaxy and that this CO component is inflowing due to 
some other reason. 

The CO component might be co nnected to the stellar oval 
found by iGarcia-Burillo et aT l (l200l inside the circumnuclear 
ring. The position angle of this oval is ~ 66°. Unlike the Ha spi- 
ral arms, its orientation matches well with the location of the CO 
component. If we draw a line through the dynamical center with 
the oval's PA, the CO component lies almost completely at the 
southern leading side of the oval (Fig. [T] bottom left). A stellar 
oval is in essence a 'fat' bar and can drive gas inward similar to 
the large-scale stellar bar on larger scales (i.e. the bars-within- 
bars scenario). The location of the CO component with respect 
to the oval is correct for this scenario, although we cannot con- 
clude inflow from our gravitational torque results. The reason 
why the influence of the stellar oval is not more prominent in 
the gas might be that its bar stre ngth and thus th e torque it ex- 
erts on the gas is low (Garcfa-Burillo et al.ll2005l ; their Fig. 8b) 
and the amount of gas reaching the influence radius of the oval 
is limited by the circumnuclear ring and its associated star for- 
mation. The CO gas we detected in the gas bridge lies in the 
region where it would feel positive torques due to the nuclear 
oval; in the averaged torque budget, as shown in Fig. [T2]this gas 
can dominate and mask the smaller amount of inflowing gas at 
these radii. Alternatively, viscosity torques may play a more sig- 
nificant role here, as was also dicussed by Garcfa-Burillo et al. 

dlool . 

This CO component could be one of the last links in the chain 
of gas inflow down to the center The large-scale bar drives most 
gas down to the ring. Based on the gravitational torques, this 
is not the smallest possible radius- only the mass accumulation 
rate becomes small. So, it is possible for gas to move further in, 
reach the sphere of influence of the stellar oval or experience 
viscosity torques, which brings it further toward the nucleus. 
The CO component we observe does not completely reach the 
dynamical center However, spatially a nd kinemati c ally, i t con- 
nects to flie central HCN detection of lKrips etaP (l2007h . The 
HCN extends over ~0.5" around the dynamical center, and our 
CO component extends to about this radius from the center The 
detected HCN also has a detected velocity range of + 70 km 
s with the positive side where the CO component is, which 
also has a velocity of 70 km s This concentration potentially 
completes the chain of gas inflow. 

7. Summary and conclusions 

We have presented high-resolution '^00(1-0) and '2CO(2-l) 
maps of the central region of NGC 695 1 using the IRAM PdBI 
and 30m telescopes. This galaxy exhibits significant indirect ev- 
idence for recent gas inflow. It has a large-scale bar, a circumnu- 
clear SB ring, and an AGN. Molecular gas, as traced by CO, is 
the dominant phase of the ISM in these regions. 

Investigation of the CO data cubes shows that the gas dis- 
tribution can be well explained by gas response to a large-scale 
stellar bar. We further quantify this scenario by reproducing the 
observations using a model where the gas is solely responding to 
a large-scale stellar bar. The model used here has the advantage 
that it provides direct predictions for the radial motions of the 
gas. Our best model is able to reproduce the main characteristics 
of the observed CO morphology and kinematics. 

Gravity torques and gas mass accumulation rates that follow 
from this model we re computed and com pared to pre vious grav- 
ity t orque maps by iGarcfa-Burillo et al.l (|2005) and iHaan et al.l 
(l2009i) . Their method is complementary to ours and we explain 



the diflrerences between the two results. As expected, we find that 
the large-scale bar eff'ectively transports gas inward to r~400pc, 
but it no longer dominates the gravitational torque budget on 
smaller radii. 

The stellar oval reported by iGarcfa-Burillo et all (l2005b is 
a likely candidate to take over gas transport inside the circum- 
nuclear ring. Detailed investigation of the data cubes revealed 
double-peaked CO profiles at several positions inside the gas 
ring. Most notably, we detect a second CO complex near the dy- 
namical center, with a velocity off'set of Av ~ 80 km s ' , which is 
not on a stable orbit. From simple geometric and physical argu- 
ments, we conclude that this cloud is in the disk of the galaxy, is 
inflowing toward the nucleus, is likely driven by the stellar oval, 
and is indeed the last step in bringing gas close to the nucleus. 
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Appendix A: Models' goodness of fit 

In Sect. 14.21 we constructed a bar model based on an analytic potential, and derived a coiTesponding model data cube. The best 
model was estimated by comparing the model's channel maps and intensity map against their observed counterparts. In Fig. lA. II we 
present overlays of the observed and modelled spectra at different positions in the cubes to further show the goodness of fit. 

Fig. lA. II shows that the model (black dashed lines) follows the observed line emission (solid lines) very well. There are some 
local differences in the intensity. These differences occur mainly at the locations of the 'twin peaks', (-1.2,5.8) and (-1.2,-5.5). There 
the observed emission is larger than we reproduce with the model. A rejected model, with slightly different values for the parameters 
char, length (210pc), char, velocity (220 km s"'), bar strength (-0.030) and bar pattern speed (60 km s"' kpc"'), is shown with grey 
dashed lines. This model was rejected because it overpredicts the intensity in the ring. 
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Fig. A.l. CO(l-O) emission spectra at different positions in the observed data cube (solid line), the best model cube (black dashed 
line) and a rejected model cube (grey dashed line). The spatial positions in the galaxy for each spectram, in arcseconds, are given in 
the subpanels. The model shows a good agreement with the observations. The positions (-1.2,5.8) and (-1.2,-5.5) correspond to the 
'twin peaks' described in Sect. 13. II 



